function x=randpdf_2D(p1,p2,px,py,dim)
% RANDPDF
%   Random numbers from a user defined distribution
%
% SYNTAX:
%   x = randpdf(p1, p2, px, py, dim)
%       randpdf(p1, p2, px, py, dim)
% 
% INPUT:
%   p1  - probability density in 1st dimension,
%   p2  - probability density in 2nd dimension,
%   px  - values for probability density in 1st dimension,
%   py  - values for probability density in 2nd dimension,
%   dim - dimension for the output matrix.
%
% OUTPUT:
%   x   - random numbers. Run function without output for some plots.
%
% DESCRIPTION:
%   x = randpdf_2D(p, px, dim) returns the matrix of random numbers from
%   probability density distribution defined in p and px. p are the density
%   (the y axis) and px are the value (the x axis) of the pdf. p and px
%   must be of the same length.
%   dim define the output matrix dimensions, for example dim=[100 3] define
%   the 100x3 two dimensional matrix with 300 random numbers.
%
%   REMEMBER: This is not a really random number generator but only
%   some kind of transformation of uniformly distributed pseudorandom 
%   numbers to desired pdf!
% 
% EXAMPLE 1:
%   Generation of normal distributed random numbers. This is not typical
%   normal distribution because is limited from the left and right side,
%   e.g. 0 < px < 80 .
%   
%   px=0:80;
%   p=1./(10*sqrt(2*pi))*exp((-(px-40).^2)./(2*10^2));
%   randpdf(p,px,[10000,1])
%
%
% EXAMPLE 2:
%   Generation using user defined pdf.
%   
%   px=[1 2 3 4 5 6 7 8 9];
%   p= [0 1 3 0 0 4 5 4 0];
%   randpdf(p,px,[50000,1])

% By Adam Nies�ony, Opole University of Technology, Poland

% check the number of input
error(nargchk(5, 5, nargin))

% vectorization and normalization of the input pdf
px=px(:,:);
py=py(:,:);
p1=p1(:,:)./trapz(px,p1);
p2=p2(:,:)./trapz(py,p2);

% interpolation of the input pdf for better integration
% in my opinion 10000 point is sufficient...
pxi=[linspace(min(px),max(px),10000)]';
pyi=[linspace(min(py),max(py),10001)]';
pi1=interp1(px,p1,pxi,'linear');
pi2=interp1(py,p2,pyi,'linear');

% computing the cumulative distribution function for input pdf
% assuming that the 2 dimensions are independent, the 2D cumulative
% distribution is equal to the product of the independent cumulative
% distribution functions.
cdfp1 = cumtrapz(pxi,pi1);
cdfp2 = cumtrapz(pyi,pi2);


% finding the parts of cdf parallel to the X axis 
ind1=[true; not(diff(cdfp1)==0)];
ind2=[true; not(diff(cdfp2)==0)];


% and cut out the parts
cdfp1=cdfp1(ind1);
cdfp2=cdfp2(ind2);
pi1=pi1(ind1);
pxi=pxi(ind1);
pi2=pi2(ind2);
pyi=pyi(ind2);

% generating the uniform distributed random numbers
uniformDistNum1=rand(dim);
uniformDistNum2=rand(dim);

% and distributing the numbers using cdf from input pdf
userDistNum1=interp1(cdfp1,pxi,uniformDistNum1(:)','linear');
userDistNum2=interp1(cdfp2,pyi,uniformDistNum2(:)','linear');

% making graphs if no output exists
if nargout==0
    subplot(3,4,[1 2 5 6])
    [n,xout]=hist(userDistNum,50);
    n=n./sum(n)./(xout(2)-xout(1));
    bar(xout,n)
    hold on
    plot(pxi, pi1./trapz(pxi,pi),'r')
    hold off
    legend('pdf from generated numbers','input pdf')

    subplot(3,4,[3 4 7 8])
    plot(pxi, cdfp,'g')
    ylim([0 1])
    legend('cdf from input pdf')

    subplot(3,4,[9:12])
    plot(userDistNum)
    legend('generated numbers')
else
    x(:,1)=reshape(userDistNum1,dim);
    x(:,2)=reshape(userDistNum2,dim);
end